function [corr0_dX_m_n_cond_on_z,Cov0_dX_m_n_cond_on_z] = CalcTimeSeriesCorr1_cond_on_z(X_Nxz,prob_Nxz,Phi_Nxz,m,n,...
    Phi_Nxz_m,Phi_Nxz_n)
% computes 
% corr(X(t+m)-X(t),X(t+m+n)-X(t+m)|z(t)=z) = A/(sqrt(B)*sqrt(C))
% A = cov(X(t+m)-X(t),X(t+m+n)-X(t+m)|z(t)=z)
%   = cov(X(t+m),X(t+m+n)|z(t)=z)-var(X(t+m)|z(t))
%     - cov(X(t),X(t+m+n)|z(t)) + cov(X(t),X(t+m)|z(t))
% B = var(X(t+m)-X(t)|z(t)=z)
%   = var(X(t+m)|z(t))+var(X(t)|z(t))-2*cov(X(t+m),X(t)|z(t))
% C = var(X(t+m+n)-X(t+m)|z(t)=z)
%   = var(X(t+m+n)|z(t))+var(X(t+m)|z(t))-2*cov(X(t+m+n),X(t+m)|z(t))
    
    [NN Nx Nz] = size(prob_Nxz);

    if nargin<=5
        Phi_Nxz_m = mpower2(Phi_Nxz,m);
    end
    
    if nargin<=6
        Phi_Nxz_n = mpower2(Phi_Nxz,n);
    end
    
    Phi_Nxz_mn = Phi_Nxz_m*Phi_Nxz_n;
    
    prob_z = reshape(squeeze(sum(prob_Nxz,[1 2])),[],1);
    
    % E[X(t)|z(t)] and Var[X(t)|z(t)]
    E0_X0 = reshape(squeeze(sum(prob_Nxz.*X_Nxz,[1 2])),[],1)./prob_z;
    Var0_X0 = reshape(squeeze(sum(prob_Nxz.*(X_Nxz.^2),[1 2])),[],1)./prob_z - E0_X0.^2;
    
    % E[X(t+m)|z(t)] and Var[X(t+m)|z(t)]
    E0_Xm = reshape(squeeze(sum(prob_Nxz.*reshape(Phi_Nxz_m*X_Nxz(:),[NN Nx Nz]),[1 2])),[],1)./prob_z;
    Var0_Xm = reshape(squeeze(sum(prob_Nxz.*reshape(Phi_Nxz_m*(X_Nxz(:).^2),[NN Nx Nz]),[1 2])),[],1)./prob_z - E0_Xm.^2;
    
    % E[X(t+m+n)|z(t)] and Var[X(t+m+n)|z(t)]
    E0_Xmn = reshape(squeeze(sum(prob_Nxz.*reshape(Phi_Nxz_mn*X_Nxz(:),[NN Nx Nz]),[1 2])),[],1)./prob_z;
    Var0_Xmn = reshape(squeeze(sum(prob_Nxz.*reshape(Phi_Nxz_mn*(X_Nxz(:).^2),[NN Nx Nz]),[1 2])),[],1)./prob_z - E0_Xmn.^2;
    
    % Cov(X(t),X(t+m)|z(t))
    Cov0_X0_Xm = reshape(squeeze(sum(prob_Nxz.*X_Nxz.*reshape(Phi_Nxz_m*X_Nxz(:),[NN Nx Nz]),[1 2])),[],1)./prob_z - E0_X0.*E0_Xm;
    
    % Cov(X(t),X(t+m+n)|z(t))
    Cov0_X0_Xmn = reshape(squeeze(sum(prob_Nxz.*X_Nxz.*reshape(Phi_Nxz_mn*X_Nxz(:),[NN Nx Nz]),[1 2])),[],1)./prob_z - E0_X0.*E0_Xmn;
    
    % Var[X(t+m)-X(t)|z(t)]
    Var0_diff_Xm_X0 = Var0_Xm + Var0_X0 - 2*Cov0_X0_Xm;
    
    % Cov(X(t+m),X(t+m+n)|z(t))
    Cov0_Xm_Xmn = reshape(squeeze(sum(prob_Nxz.*...
        reshape(Phi_Nxz_m*(X_Nxz(:).*(Phi_Nxz_n*X_Nxz(:))),[NN Nx Nz]),[1 2])),[],1)./prob_z - E0_Xm.*E0_Xmn;
    
    % Var[X(t+m+n)-X(t+m)|z(t)]
    Var0_diff_Xmn_Xm = Var0_Xmn + Var0_Xm- 2*Cov0_Xm_Xmn;
    
    Cov0_dX_m_n_cond_on_z = Cov0_Xm_Xmn - Var0_Xm - Cov0_X0_Xmn + Cov0_X0_Xm;

    corr0_dX_m_n_cond_on_z = Cov0_dX_m_n_cond_on_z./(sqrt(Var0_diff_Xm_X0).*sqrt(Var0_diff_Xmn_Xm));

end